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Abstract. We consider the fluctuations of generalized currents in stochastic 
Markovian dynamics. The large deviations of current fluctuations are shown to obey a 
Gallavotti-Cohen (GC) type symmetry in systems with a flnite state space. However, 
this symmetry is not guaranteed to hold in systems with an inflnite state space. A 
simple example of such a case is the Zero- Range Process (ZRP). Here we discuss in 
more detail the already reported [1] breakdown of the GC symmetry in the context 
of the ZRP with open boundaries and we give a physical interpretation of the phases 
that appear. Furthermore, the earlier analytical results for the single-site case are 
extended to cover multiple-site systems. We also use our exact results to test an 
efflcient numerical algorithm of Giardina, Kurchan and Peliti j2] , which was developed 
to measure the current large deviation function directly. We flnd that this method 
breaks down in some phases which we associate with the gapless spectrum of an effective 
Hamiltonian. 
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1. Introduction 

An important step in the understanding of non-equilibrium systems has been the 
development of so-called fluctuation theorems [HI [4] which quantify irreversibility by 
relating the probability of some "backward" process to that of a corresponding "forward" 
process. Formally, these fluctuation relationships all derive from a statement about the 
relative probabilities of a trajectory in phase space and the time-reversed trajectory. 
The various theorems found in the literature can be divided into two broad classes: 
transient relations (which are exact for finite times) and steady-state relations (which 
hold only asymptotically in the long-time limit) — see, e.g., [5], [6] for recent discussion 
of the various interconnections. In the present paper we will mainly be interested in 
asymptotic symmetries, specifically whether (or not) a given quantity satisfies a relation 
of the form 

^(-'■•*' - e-« (1) 



Here V{r,t) is the probability to observe, over time interval [0,t], a time-averaged value 
r = R/t of some time-extensive quantity R (e.g., entropy production). The symbol 
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~ means logarithmic equality in the limit of large t. In this paper, we refer to the 
relationship ^ as "the Gallavotti-Cohen fluctuation theorem" or simply "the fluctuation 
theorem". 

Historically, the concept of a fluctuation theorem emerged from computer 
simulations investigating the entropy production fluctuations in a sheared fluid |7j. A 
rigorous derivation of the form ([T]) was then given for the steady state of deterministic 
systems [HI E] where the entropy-like functional r can be identified with the rate of 
phase space contraction. Later proofs addressed stochastic dynamics [IOl[II], leading to 
a general property of the large deviation function sometimes referred to as "Gallavotti- 
Cohen symmetry" (or GC symmetry for short). Here r is associated with a current 
of some quantity through the system (e.g., the particle current in a lattice gas) and, 
for bounded state-space, the relation ([T]) holds for arbitrary initial condition. Recent 
work by Kurchan [l2] has also explored how to recast the deterministic fluctuation 
theorem as the vanishing-noise limit of the stochastic one. In parallel to the theoretical 
development there has been much successful work on experimental verification of 
fluctuation theorems; for reviews see [U [131 El- In particular, we note that recent 
experiments using an optical defect in diamond provide a simple realization of a two- 
state stochastic system [T5l [T6]. 

Very recently there has been considerable theoretical, experimental and numerical 
interest in cases where the symmetry ([T]) breaks down, see e.g., [TT l fTS l fT9l[20l [Tl [21 1 [22] . 
It is now understood that this effect can be attributed to certain boundary terms which 
become relevant in the case of infinite state space; see section [2] below for a detailed 
exposition in the context of stochastic Markovian dynamics. For deterministic systems, 
the effect of unbounded potentials was discussed by Bonetto et al. [20]. They argued for 
the restoration of the symmetry by removal of the "unphysical" singular terms. Earlier 
a similar phenomenon was found for a model of a trapped-Brownian particle treated 
via a Langevin equation [I7l [23]. Within the Langevin framework, discussion of some 
related subtleties can also be found in [241 [25] . 

In [I] we provided a general discussion of this breakdown of Gallavotti-Cohen 
symmetry for many-particle stochastic dynamics within the context of a particular 
model — the zero-range process (ZRP). The ZRP plays a paradigmatic role in non- 
equilibrium statistical mechanics, hence such work offers a contribution to general 
understanding as well as providing concrete results for an important model. In the 
present paper, we give details of the analytical calculations leading to those results 
and compare them with direct evaluation of the current large deviation function using 
the recent algorithm of Giardina, Kurchan and Peliti [2]. We also discuss some 
generalizations and present new results for the multi-site ZRP. 

Interestingly, in the zero-range process, we do not find a constant value for the ratio 
of probabilities for large forward and backward current fiuctuations. This is in stark 
contrast to analytical arguments [TTl l20l l25| and numerical work [21] for other models. 
We argue below that the failure to observe this form of "extended fiuctuation theorem" is 
due to strong correlations in our model between the boundary terms and the integrated 
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current. These correlations persist even in the long-time limit; it would be interesting 
to see if effects of this type can be observed in any experimental situations. 

The plan of the rest of the paper is as follows. First, in section ^ we give a 
general derivation of the Gallavotti-Cohen fluctuation theorem for stochastic systems 
and indicate also the relation to transient fluctuation theorems. Then, in section [3l we 
introduce our zero-range model and outline its treatment within this general framework, 
section [4] is devoted to a detailed calculation of the current large deviation function for 
the single-site version of this model, giving a concrete example of the breakdown of 
GC symmetry. In section [5] the analytical approach is complemented by use of the 
algorithm of [2] to obtain new numerical results for the large deviation function (and 
general insight into the applicability of the "cloning" method). Significantly, in section [6] 
we extend our discussion to larger systems, demonstrating how information about the 
current fluctuations in an L-site ZRP can be obtained by mapping to an effective single- 
site model. Section [7] contains some conclusions and general perspectives. 

2. The fluctuation theorem for Markovian dynamics 

2.1. Central argument 

Here we present a derivation of the fluctuation theorem for generalized currents of 
Markov processes defined on a finite configuration space. Our argument is based on 
that of [n]. 

Consider a continuous time Markov process which satisfies detailed balance in 
the stationary state. The system can be described by the transition rates w'^J}^ from 
configuration a to a'. In addition, consider a counter J, the value of which increases 
by the amount 9^/^ at each transition a a'. Here the matrix 0, which is required 
to be real and antisymmetric, can describe any type of real or abstract current in the 
system. As an example one can consider the particle current through a given bond: in 
this case Q^'a is the number of particles hopping across the bond at a transition cr — > a' 
(which can be positive or negative depending on the direction of hopping). At t = 
the counter J is set to zero; for any positive time it is a functional acting on the paths 
(sequences of configurations) from time to time t and can be written as 



Here n is the number of configurations visited during time t. In the following we 
refer to J as the "time-integrated current" (or, where no confusion can arise, simply as 
the "current"). Note that, due to detailed balance, the mean of this current {Jit)) is 
always zero in equilibrium. 

We define the driven system by the modified transition rates 



j{t, {a}) = E e, 



(2) 



k=l 



(3) 
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from configuration a to a', where -E" is a driving field conjugated to the specific current 
under consideration. In what follows we show that in this driven system the probability 
distribution function V{J,t) of the random variable J{t) satisfies the relation 



V{-J,t) 

asymptotically for large times, provided the state space (i.e., the number of possible 
configurations) is finite. Here the power EJ can be interpreted as the work done on the 
system by the external field. This is the statement of the fiuctuation theorem. 

In cases where the sum J2n 0cr„,CT„+i gives zero for all periodic paths cr„ in the 
configuration space, J becomes a simple function of the initial and final configuration, 
i.e., there is no real dependence on the history. In this special case not only the 
original but also the above-defined "driven" system would satisfy detailed balance in 
the stationary state. In the following we assume that this is not the case, i.e., there are 
periodic paths in the configuration space for which the sum J2n^(Tn,(Tn+i is non-zero. 

For a given non- equilibrium model with rates Wa'a one can apply a reversed 
argument. In this case physically one can think of -E as a negative driving field, 
conjugated to the current under consideration, which is needed in order to "restore" 
detailed balance. If, for a specific value of E, the rates w'^^ (defined by (l3|)) lead to 
detailed balance then the fluctuation relation (Ilj) holds for this current. This gives some 
freedom for the quantity J which enters the fluctuation relation. The action functional 
of [H] {W in that paper) corresponds to the specific choice Q^'a = In with E = 1, 
which leads to w^?^ = a/uvv^W^- These rates indeed lead to detailed balance, since for 
each pair of configurations the forward and backward transition rates are equal. This 
also implies that in the corresponding equilibrium distribution every configuration has 
the same weight. 

Let us now define the rate w({cr}) for a full path as 

_tO_ *-*n-l 

^(l^l) = ^<xiaoW<72<7i ■ ■ -u^ana^-ie "oe n ...e , (5) 

where tk denotes the time when the transition from configuration ak to cxfc+i happened 
and Tk = (Z]<j' ^CT'cTfe) ^ corresponds to the overall exit rate from configuration cxfc. The 
conditional probability of such a path with transition times between tk and tk + dt, 
provided at time the system starts in configuration ctq is w{{a})dt"-~^ . Using (l3|) one 
can readily show that 



pZw ({a}-) ' 

where {aY'^^ is the time-reversed path of {a} and is the equilibrium probability of 
configuration a. In reference |26| the above relation is referred to as the "non-equilibrium 
detailed balance condition" with EJ being the work done on the system. This also 
suggests the identification of EJ as the work. 

We note that in the case of discrete time dynamics the above scenario is very 
similar, the only difference is that here the quantities Wfj'a and w{{a}) denote transition 
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probabilities rather than rates and the exponential factors in ([5]) are replaced by diagonal 
elements of the transition matrix w. 



2.2. Proof of the asymptotic fluctuation theorem for finite systems 

In the calculation we use the so-called quantum Hamiltonian formalism where a basis 
vector is associated with each possible configuration and the state of the system (a 
probability measure on the configuration space) is denoted by a vector \P) in this space 
with the normalization {s\P) = 1. Here (s| is a row vector with components (1, 1, 1, . . .). 
In this formalism the master equation takes the form 

||P) = -H\P), (7) 

which is similar to the Schrodinger equation in imaginary time. The transition rates are 
in the off-diagonal elements of H and the conservation of probability requires 

{s\H = 0. (8) 

For more details on this formalism see \27\ . 

As a first step of the proof we introduce the joint probability distribution function 
Va{J, t), which denotes the probability of being in configuration a and having the value 
J of the current at time t. In what follows we calculate the generating function (e"'^"'^*^). 

{e-''^'^) = {s\g{t)) , (9) 

where 

g{t)^ = Y.'PMt>~''- (10) 
J 

The time derivative of g{t) is 

= EE^-'(^''^)^-'e"'''e-'''--' -^7(t).E^-v 

J' cr' a' 
a' _ a' 

= ~H{\)„„,g{t),., (11) 

where H{X) is a modified Hamiltonian in which the transition rates corresponding to 
a' ^ a are multiplied by the factor e'^'^""' . Note that H{X) is a non-stochastic matrix 
with {s\H{\) 7^ (s| for A 7^ 0. Since |5'(0)) is identical to the initial measure \Pq) the 
generating function takes the form 



-\j{t)\ 



Po). (12) 



The long-time behaviour of this quantity is characterized by the function 

e(A) = -^lim iln(e-^-'W). (13) 
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One finds from (fT2ll that as long as the configuration space is finite, e(A) is given by the 
lowest eigenvalue eo(A) of i/(A), since 

(e-^^W) = Y: mPo) e-^-^^^* ^ (^l^o) (HPo) e-'"^'^', (14) 

i 

Where ej(A) are the eigenvalues and ipi are the eigenvectors of if(A). It is easy to show 
that e(A) is related to the large deviation function 

e{j) = -\imhnV{tj,t) (15) 

t^co t 

of the time-averaged current j = J/t through a Legendre transformation: 

e(j) = max{e(A) - Xj}, e(A) = mm{e(j) + Aj}. (16) 

As a second step in the proof we show that H has the symmetry property 

HiXf = P-^'H{E - A)Peq, (17) 

where Peq is a matrix with the equilibrium probabilities on the diagonal and zero 
elsewhere. For the rhs of (fTTI) one finds 



p-'HiE-\)P,^ , =_^^,^e-(^-^)®'^''^^^(l-(5.v) + E^P-'^-'-' (18) 



Po 



p 



which leads to 



^..'e-^®^^'=^.ve-(^-^)®^'^%^ (19) 



for the non-diagonal elements of equality (jlTl) . Using ([3]) and the fact that the matrix 
G is antisymmetric, the above condition takes the form 

^ll'V^ = Klf^ (20) 

which is just the detailed balance condition for the equilibrium system and is trivially 
satisfied. This proves the relation (fTTI) (The diagonal part of both the Ihs and rhs is 
just the diagonal part of the original non-equilibrium Hamiltonian.) 

As a corollary of the symmetry property (fTTll one finds the important relation 

e(A) = e(E-A). (21) 
For the large deviation function e(j), using the Legendre transformation ( |T6l ) this implies 

ei-j)=Ej+eij), (22) 
which is equivalent to ([Ij). 
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2.3. Transient fluctuation theorem 

The relation (JH) is true only asymptotically, for large times. However, due to the 
symmetry (ITTl ). one can show that for specific initial conditions it can be made exact for 
any finite time. Namely, taking the equilibrium distribution {p^"^) as the initial condition 
then for the generating function one finds 

= (s|e-^(^-^)>^<^) = (e-(^-^)-^W). (23) 

This implies that the relation ([4]) holds exactly for any finite time (but only for this 
specific initial condition). This is the statement of the transient fluctuation theorem of 
Evans and Searles [28l[29l[3]. 

2.4. Boundary terms 

As mentioned above, in a non-equilibrium system one has some freedom to choose the 
current to be considered. Here we investigate whether two different choices really give 
two independent relations for the fluctuations of the non-equilibrium system. In order 
to satisfy the conditions of the theorem, by applying —E field in the non-equilibrium 
model one should get back to detailed balance, which can be formulated as 

"^^^^^ = (24) 

Here V„ is the energy of configuration a in the equilibrium system. Taking the logarithm 
of (l24l) and summing up for a path results in 

j*(t)-i5;j(t) = K..-Kfl„, (25) 

where J* is the action functional of [TT] and V^i^i^gj^j is the potential in the initial (final) 
state. This shows that any two currents (that satisfy the conditions of the fluctuation 
theorem) differ only in boundary terms. In finite systems these terms are bounded, 
consequently their contribution to the net current is negligible in the limit, where t ^ 00. 
In these systems the fluctuation theorems for different currents give essentially the same 
information. 

In order to demonstrate this, consider a driven exclusion process on a finite lattice. 
Here the number of possible configurations is clearly finite. There are many ways of 
defining the integrated particle current in the system: one can consider the current 
through a given bond, take a space-averaged current or measure the distance travelled 
by a tagged particle (in periodic systems). All these definitions give a different current 
for finite times but in the t (yo limit they are essentially the same, since the difference 
between them is bounded while they are proportional to t. 

However, for systems with infinite configuration space this argument does not 
hold and the boundary terms can be relevant even for large times, leading to several 
independent relations for the fiuctuations. 
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The fact that the difference between two currents can be written as a boundary 
term suggests that the same holds for the difference between two initial states. Suppose 
that the transient fluctuation theorem holds for J^^\ E^^^ with the initial measure 
eq(i) _ ^-vy\ gij^pg currents E(i)J(i) = E^^) j(2) + v;^^) - Vj^) - kJ^) + I/j^) , 

J. u yj t-'ini ^fin "-'mi '-'fin ' 

the transient fluctuation theorem holds for any current J with a boundary term as 

J,o„ = + K., - K«„ -VZ'-yZ, (26) 
where V is the equilibrium potential corresponding to the current J and V^°' = Inp™' 
with p™^ being the initial measure (Jcorr satisfies the theorem with E = \). 

2. 5. Breakdown of the fluctuation theorem in infinite systems 

For infinite systems the above argument breaks down at (fT4ll . where one identifies 
e(A) defined in (fT3l) with the lowest eigenvalue of H{\). In cases where H is infinite- 
dimensional, there is no guarantee that the scalar products appearing in ( |T4l) are finite. 
An explicit example of this breakdown is given in [Tj and will be discussed in detail 
below. Note however, that the derivation of the transient fiuctuation theorem still 
holds even for infinite systems. Consequently, a violation of the asymptotic fiuctuation 
theorem requires the unbounded growth of the boundary term in fl26l) . 

To avoid possible confusion, we remark here that in the original dynamical systems 
formulation [H [9] of the GC symmetry, a critical value r* already appears above which 
the relation (P) does not hold. This is a direct consequence of the fact that the large 
deviation function (equivalent to e in our notation) diverges outside a finite interval. 
In contrast, in the stochastic framework of this paper, we consider situations where the 
large deviation function is always finite (i.e., formally r* = oo) but the symmetry may 
still have a restricted range of validity due to the relevance of boundary terms. Indeed, 
it has also recently been understood in the dynamical systems context, that adding 
singular boundary terms to a system with finite r* leads to a further reduction in the 
symmetry regime, see [20] . 

2. 6. Extended fluctuation theorem 

The argument we present in this section is largely based on that of van Zon and Cohen 
[TTl [23] . where they describe a generic scenario for the breakdown of the fiuctuation 
theorem in the context of Langevin dynamics. This leads to an "extended fiuctuation 
theorem", where the quantity e(j) — e(— j) is linear for small j (as suggested by the 
fiuctuation theorem) but, after an intermediate crossover regime, saturates to a constant 
value for large j. Similar arguments are also given in [25]. 

Our starting point is equation (l26ll which we now write for the time-intensive 
quantities: 

Jcorr =J+t"'(Ani-5fi„). (27) 

Here jcorr is the corrected current for which the fiuctuation theorem holds, and 5 is a 
boundary term, which depends only on the initial/final configuration of the history. 
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First we assume that the probability distribution Vi^b) of B has an exponential tail 
for very large/small values of B: 

Here we set to infinity if B is bounded from above/below. Focusing on the 

case where the initial state is the steady state, one can assume that the probability 
distribution of B is identical for the initial and final state, moreover, for large 
measurement times they become independent. This implies that the probability 
distribution of 

6 = ri(5i„i-5fin) (29) 

is of the form 

P(,)(6,^)~e-*^l^ (30) 

where A = min(A+, A^). 

As a next step we assume that b and jcorr as random variables are independent. In 
this case the large deviation function of the fluctuations of j = jcorr + b are determined 

by 

e(j) = min(ec(jcorr) + ^Ijcorr - jD- (31) 

Jcorr 

Here Cc, which describes the large deviations of jcorr, satisfies the fluctuation theorem, 
i.e., ec(— j) — ec(j) = j. In figure [T], one can see how e(j) can be obtained from ec(j) 
graphically by using flSTl) . The points —ji and j2 are defined as those values of j where 
the derivative of the convex function edJ) becomes —A and A respectively. Outside the 
interval (—^1,^2) the function e(j) is linear with slope ±A. Consequently, the quantity 
e(— j)— e(j) is linear from zero to j = min(ji, ^2), and after a crossover regime it saturates 
at j = max(ji,j2) to a constant value C. We remark that in the case where edj) is 
symmetric with respect to (j) (e.g. in the case of Gaussian fiuctuations) C = 2A{j). 

In a more general situation, where the initial state is not stationary, the fiuctuations 
of B are different in the initial and final state (we still require that the tails are 
exponential). 



In this case the large deviations of b take the following form: 

g-tAi6 ^ ^ Q 



.M, . Z ^ (33) 



where Ai = mm{A^^, A^) and A2 = m.m{A^^, A'^) . This leads to a variant of the 
extended fiuctuation theorem, where the the quantity e(— j) — e(j) does not saturate 
but becomes linear with a non-zero slope for large j. 

We stress that the above argument relies strongly on the assumption that the 
boundary terms are independent of the bulk contribution. There are various examples in 
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-j2 -Jl (j) ji J2 

3 



Figure 1. A schematic plot showing the large deviation functions e(j) and ec(j). 
While the fluctuation theorem holds for ec(j), a modifled form is found for e(j), as 
plotted in the inset 

the literature where this holds, and the above extended form of the fluctuation theorem 
was found, e.g., [23l [13, El]. In contrast, in the remainder of this paper, we discuss 
a simple stochastic Markovian system where there is a strong correlation between the 
boundary and bulk terms and hence this extended fluctuation theorem is not expected 
to hold. 

3. Zero-range process 

Here we demonstrate how the general formalism of the preceding section is applied to a 
specific model — the zero-range process (ZRP). First introduced by Spitzer in 1970 [30] . 
the ZRP now plays a paradigmatic role in non-equilibrium statistical mechanics, see [3l] 
for a recent review. In particular, for certain choices of parameters the model exhibits 
a condensation transition [32l [33] in which a macroscopic proportion of particles pile 
up on a single site. Condensation phenomena are well-known in colloidal and granular 
systems [34] and also occur in a variety of other physical and non-physical contexts [31] . 

3.1. Model 

We study the one-dimensional partially asymmetric zero-range process with open 
boundaries [35j — see figure [3 Each lattice site can be occupied by any integer number 
of particles, the uppermost of which hops randomly to a nearest neighbour site after 
an exponentially distributed waiting time. In the bulk particles move to the right (left) 
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qWn PWn 



JWn J)Wn 




12 3 4 5 L-1 L 

Figure 2. Schematic representation of the ZRP on an open L-site lattice 



with rate pWn {qWn) where w„ is a function of the number of particles n on the departure 
site {wq = by definition). Note that Wn = n would correspond to free particles whereas 
all other forms represent an attractive or repulsive inter-particle interaction. The top 
particle on the leftmost lattice site (site 1) leaves the system with rate 7W„ whereas 
new particles are injected with rate a. Correspondingly, on the rightmost site (site L) 
particles are removed (injected) with rates (S). For later convenience we label each 
bond by the site at its left-hand end, i.e., the /th bond is between sites / and / + 1. 

In the quantum Hamiltonian formalism this dynamics is encoded in the Hamiltonian 



H 



1=1 



ai a;+ 1 - di) + q^afai^^ - di+i) 
+a{at - 1) + 7(ar - di) + 5{at - 1) + Pia} 



(34) 



/ wi 

W2 

W3 





where a"*" and a are infinite-dimensional particle creation and annihilation matrices 

/ 
10 
10 
10 



(35) 



and d is a diagonal matrix with the (i, j)th element given by Wi6ij. 



3.2. Stationary state, current fluctuations 

The steady state of the ZRP is given by a product measure (in other words, the 
stationary distributions for each site are uncorrelated): 

|p*) = |p*) ® \p;) ®...®\pi) (36) 

where |P;*) is the probability vector with components 

n 

P*{ni = n) = ^l[w-\ (37) 

Zj] ■ -I 
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Here the empty product n = is defined equal to 1 and Zi is the local analogue of the 
grand-canonical partition function 



n 



Zi^Z{z,) = J2zrY[^-\ (38) 

n=0 i=l 

The fugacities zi are uniquely determined by the hopping rates a, /5, 7, 5, p and q. 
However, for Wn bounded (i.e., lim^^ooWn = a with a < 00), Zi has a finite radius of 
convergence. For parameters leading to fugacities outside this radius of convergence, a 
growing boundary condensate occurs |35| . 

We are interested in the application of the fiuctuation theorem to this model, 
for the case where the boundary parameters are chosen to give a well-defined steady 
state, i.e., without boundary condensation. As discussed in the preceding section, one 
can consider a variety of different "currents" through the system. However, since the 
state space is unbounded (each lattice site can contain an arbitrarily large number 
of particles) we anticipate the possibility of relevant boundary terms and different 
fiuctuation relationships. 

A natural choice is to look at the physical current of particles across a particular 
bond. For example, suppose we choose to focus on the particle current into the system 
(i.e., across the 0th bond), then the modified Hamiltonian is given by 

^ 1=1 

+a(a+e"^ - 1) + 7(«re^ " ^^i) + ^((^t - 1) + /^(% - ^^l)}- (39) 

In the next section we will examine closely the spectrum of this Hamiltonian for the 
L = 1 single site case. Here we recap some known results for the general L-site 
obtained in [36| . 

If Wn is unbounded, (i.e., lim^^poWn = 00) then H{\) has a gapped spectrum for 



all A with lowest eigenvalue given b; 

(p-g)(e^ - 1) 

eo(A) 



a/3 (2)^"' e-^- 75 



lip-q- P) + Pip-q + i) 



(40) 



The right eigenvector l^po) corresponding to (l40l) has the same form as the stationary 
product measure ( l36ll with fugacities 

[(ae-^ + 5){p-q)- a(3e-' + 7<5] (f)'^' - 7^ + a/^e"^ (f)^'' ^ ^ 

zi = — — . (41) 

l{p-q-f3) + f3{p-q + 7) (1) 

and the left-hand eigenvector is also a product state with one-site marginal "fugacities" 

_ /?7(e^ - 1) (f )^"' + ie\p - g - /3) + /5(p - g + 7) (h)"^"' 

— • (42) 

l{p-q-f3)+f3{p-q + 7) (2) 

I The leading order term in an L ^ 00 expansion for the bulk-symmetric case, corresponding to the 
limit p = q = 1 in (|40l) . has also been obtained by an additivity principle [37] and by field-theoretic 
methods [38]. 
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Figure 3. Schematic representation of the single-site ZRP with Wn = \. 



Note that this result is independent of the details of Wn and one can also show that 
it is the same for currents across all bonds meaning that the current fluctuations 
are spatially homogeneous in the long-time limit. Furthermore the scalar products 
appearing in equation f |T4ll are flnite and the standard Gallavotti-Cohen fluctuation 
theorem is recovered. Physically, unbounded Wn means that there is no chance for the 
temporary accumulation of a large numbers of particles on any site. 

However, ij Wn is bounded, then there is a crossover to a gapless spectrum at 
some value of A and also the scalar products in equation (fT4l) can be infinite. In fact, 
{iPq\Po) and {s\ipo) are related to the distribution of initial and final boundary terms 
and the condition for (ipolPo) to diverge obviously depends on the initial state. Thus, 
although the symmetry of the eigenvalues remains, we expect the breakdown of the usual 
Gallavotti-Cohen fluctuation theorem in a spatially homogeneous and initial-condition- 
dependent way. To explore this issue in more detail we present, in the next section, 
explicit calculations for the single-site case. 

4. Analytical results for single-site model 

As a simple example of the ZRP with bounded Wn, we now focus on the single-site model 
with Wn = 1 — see figure [3l Even this apparently simple model exhibits a rich phase 
behaviour, as shown in [1]. One must now distinguish between the particle currents 
across only two bonds viz. the 0th ("input") and 1st ("output"). Since fiuctuations across 
the two bonds are simply related by left-right refiection, we will consider only the former 
except where explicitly stated otherwise. Note that the occupation number n of the site 
performs a random walk on a semi-infinite lattice but with two independent processes 
for movement in each direction. For a well-defined steady state this random walk should 
be biased towards the reflecting boundary at n = 0, i.e., we require a + 6 < (3 + •y. If 
this condition is not met then there is a growing condensate on the site. 
Let us first consider an initial Boltzmann distribution given by 

oo 

\Po) = {l-x)Y.x^\n) (43) 

ra=0 

where \n) denotes the configuration with site occupied by n particles (i.e., a column 
vector with a T in the nth position and 'O's elsewhere) and the "fugacity" x is less than 
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Values of A 



/3+7-5 

27<5 



_ <5-/3a;2 + ^(5-/3x2)2+4Q:7x2 



gA4 ^ /3(l-^)+7 
7a; 



27x2 



Corresponding values of j 

■ _ (/3+7-^)2-a7 . _ /3(/3+7-(5)2-a75 
Ja — 5 — (/3+^)(/3+7-5) 

— /3+7 



-((5-/3a:2) 

X 

■ ^ a/37x2-5[/3(l-x)+7]^ . ^ a7-[/3(l-x)+7]^ 

"'e — a::(/3+7)[/3(l-x)+7] ' -If " /3(l-x)+7 



Table 1. Transition values for input current fluctuations in single-site ZRP. 



1 for normalization. For example, the choice x = (a + + 7) is the steady-state 

initial condition. Note also that the limit x — > corresponds to the empty-site case and 
that by ergodicity this gives the same result as any fixed particle configuration. 

To obtain the large deviations of input current, we need to calculate the matrix 
element (s|e~^*|Po) where the modified Hamiltonian H{X) is given by 

/ a + 5 -76^ -p . . . \ 

-ae-^ -6 a + /3 + 7 + 5 -76^ - /3 

-ae~^ - S a + P + 'j + 5 -76^ - /? ... (44) 

-ae"^ -6 a + (3 + + 6 ... 

V ^ '■ J 

In appendix lA.ll we first present an explicit calculation of the spectrum of H{X) and 
thereby illustrate the crossover to a gapless regime. Then, in appendix IA.21 we show 
how to obtain an integral representation of the matrix element (s|e^*^^''|Po) and thus 
extract the large deviation behaviour. This enables us to construct and explain the 
"phase diagram" for current fluctuations as shown in section 14. 1[ Finally, in 14.31 we 
discuss some straightforward generalizations of these results. 



4.1. Phase diagrams 

In order to extract the large-time behaviour from the integral representation (IA.13I) we 
can use the method of steepest descents with saddle-point at z = 1. However, care 
must be taken due to the poles in the integrands (at z = (0x)^^,0 and y). If 
the saddle-point contour has to be deformed through one of these poles then we must 
take into account the contribution of the residue at that pole. For a given A, e(A) is 
then determined by the term which decays most slowly with t. This yields changes in 
behaviour at the values of A given in table [H where we have defined for convenience the 
parameter combination 

7] = ^[{(3 + 7)2 -(36- a7]2 - Aa(3-f6. (45) 

Phase A: A < Ai (see figured]). Here e(A) takes the form 

e(A) = a(l -e-^) +7(1 -e^) , (46) 
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Figure 4. The x-X phase diagram for a = j = 5 = 0.1, /? = 0.2. Red hne: A = Ai, 
green Hne: A = A2, blue hne: A = A3(a;), cyan hne: A — X4{x). 



Phase B: 



Phase C: 



Phase D: 



which does not coincide with the lowest eigenvalue of H. The reason is that 
here the quantity {s\ipo) in (fT4l) diverges. This phase corresponds to large 
forward currents. 

[(x < Xc) A (Ai < A < A2)] V [{x > Xc) A (Ai < A < A4)], where we defined Xc 
as 

_ -v + i/3 + 7)' -ai + /36 

W{P + 7) ^ ^ 

(see figured]). In this region the spectrum of H is gapped and e(A) coincides 
with the lowest eigenvalue: 



e(A) 



(48) 



/3 + 7 /5 + 7 

The quantities {s\%Ijq) and (V'ol-fo) are finite. 

(x < Xc) A (A2 < A < A3) (see figured]). In this phase the spectrum of H is 
gapless. The scalar products {s\%Ijo) and (t/'oI-Po) are finite (here \iPq)^ which 
is the ground-state of if, has to be understood as l^/'fc^o)), consequently e(A) 
is given by the lowest eigenvalue: 

e(A) = a + (3 + -i + 5- 2^{ae-^ + 5) (/? + 7e^) (49) 

[(x < Xc) A (A > A3)] V [{x > Xc) A (A > A4)] (see figured])- In this phase the 
quantity (^o|-Po) diverges and e(A) differs from the lowest eigenvalue of H{\): 

ae~^ + 6 



e(A) = a + /5 + 7 + (5-/5 + 7e 



X 



X 



(50) 



It is interesting that here the large deviation of current fiuctuations retains 
a dependence on the initial state of the system. 
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Figure 5. The x-j phase diagram for a = 7 = (5 = 0.1, /? = 0.2. At first order 
transition fines (A-B and B-D) mixed phases appear. Red fine: j = ja or j = jb, 



green fine: j 



3c 



blue fine: 



3 



3d{x), cyan fine: j = je{x) or j = J/(a;). The 



horizontal dotted line shows the mean current, the vertical dashed line indicates the 
the specific value of x which corresponds to the steady state initial condition. 



For the physical interpretation of these phases it is better to consider e(j), which 
can be obtained from e(A) by a Legendre transformation. e(j) has the following forms 
in the different regions of figure [5l 



e(j) 



/i(«,7) 

Jj V/3+7' /9+7/ 

f,{a,^)+P{l-x) + 5{l-x- 



+ j In X 



A 
B 

C 
D 



(51) 



with 



fj{a, h) = a + b — Jp + Aab + j In 



j + y/f + Aab 
2a ■ 



(52) 



We remark that fj{a,b) has a simple physical meaning. Consider a single particle 
performing a simple random walk on the infinite one-dimensional lattice with right 
(left) hopping rate a (6) . Then the large deviation function of the distance travelled by 
this particle is given by fj{a,b). 

In phase A only the rates for the first bond determine the large deviation function. 
For such large currents particles typically pile up on the site, which then acts as an 
infinite reservoir. In this case the behaviour of the two bonds decouples, so the current 
distribution across the input bond is entirely controlled by the two Poisson processes 
at rate a and 7. This explains the first line in ([5T]) . The fact that in these (very 
unlikely) realizations the occupation number increases with time, is consistent with the 
divergence of {s\ipo)- In stochastic systems such divergence is usually associated with 
the lack of a steady state. Here, the model does have a steady state so that realizations 
involving piling-up of particles (condensation) are never observed in the infinite time 
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limit. However, as indicated, they characterize the fluctuations that can be observed 
over a large but finite measuring period t. We therefore use the term "instantanous 
condensate". 

In phase B, which always contains the mean steady-state current (A = 0), the four 
rates enter symmetrically in the large deviation function. Note that the combinations 
a/9/(/3 + 7) and 75/(/3 + 7) are the effective renormalized hopping rates of two exclusion 
particles with left (right) hopping rates a (7) and /3 (5), in the case where they form a 
bound state [39|. In this regime the occupation of the site, which maps to the distance 
between the two exclusion particles, remains finite. The existence of this "two-particle 
bound state" is also manifested in the gapped spectrum of H. 

Although the naive approach (i.e. identifying e(A) with the lowest eigenvalue of 
H) still works in phase C, the large deviation function takes a different form here. 
For these large negative currents (and small initial occupation) it is needed that the 
current across the other bond also takes a large negative value (unlike in phase A). It 
can be seen in f ISTil that the contribution of the two bonds factorize suggesting that 
they act independently. This can be traced back to the change in the spectrum which 
becomes gapless in this regime, corresponding to a two-particle unbound state. The 
distance between the two particles (or equivalently the occupation of the site) grows as 
the square-root of time in such realizations — again we refer to this temporary piling-up 
as instantaneous condensation. 

Phase D is in some sense the counterpart of phase A. Typical realizations 
contributing to these exponentially small probabilities start at a very high occupation 
(initial condensate) which then decreases during the observation time. For fixed j, this 
is possible only for sufficiently large values of x. This initial singularity is indicated by 
the divergence of {ipo\Po). 

One can see in the x-j phase diagram (see figure [5|) that in between phases A-B 
and B-D transition regions appear. These regimes correspond to a single transition line 
in the x-\ diagram along which the derivative of e(A) is discontinuous. This is entirely 
analogous to ordinary equilibrium phase transitions where different thermodynamical 
potentials are Legendre transforms of each other and in certain phase diagrams mixed 
regions (e.g. liquid-gas) appear at first order phase transitions lines. For the full analogy 
one can consider the following identifications: j (specific) volume; A ^pressure; 
e(j) Helmholtz free energy (density); e(A) — * Gibbs free energy (density). One can 
immediately see that the analogue quantity of the system size N is the measurement 
time t in our model which must diverge if a true phase transition is to exist. In the 
mixed regions of the phase diagram the system segregates in time, i.e., for some finite 
fraction of the whole measurement time the system behaves as being on one boundary 
of the mixed phase while for the rest of the time it switches to the other boundaryJ§| 

§ Mathematically speaking, knowledge of e(A) is only sufficient to determine e{j) where e(A) is 
differentiable. If e(A) is non-differentiable then e(j) is, in general, non-convex; Legendre transform 
of the non-differentiable points in e(A) then yields straight-line sections of the convex envelope of e(j). 
However in our case, the physical argument based on phase separation in time, indicates that the Hnear 
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One can formally consider t as a space dimension. Then a path in the configuration 
space of the original model becomes a configuration of the new model, see e.g., [43] . 
If the original model has finite number of configurations then this leads to only one 
(diverging) dimension (t) in the new model hence no phase transition is possible in 
this case. Therefore for a real phase transition to take place one needs infinitely many 
possible configurations in the original model. In our single site model this is achieved 
in the simplest way. 

The above analogy between ordinary equilibrium phase transitions and transitions 
in the large deviation function of current fiuctuations suggests that the ordinary 
free energy can be considered as the large deviation function of density fiuctuations, 
which is indeed the case. Vice versa, the large deviation of current fiuctuations e(j) 
can be interpreted as some kind of "dynamical free energy". The main difference 
between the two cases is that the pressure can easily be controlled and measured in 
an experiment, whereas that is not the case with A. So far it seems that the only way 
to make measurements for e(A) or e(j) is the highly inefficient method of measuring 
the probabilities of large current fiuctuations in a system where A = (zero pressure). 
However in the next section, we will see that with the help of a neat trick A becomes 
adjustable in computer simulations. Still, the question of the (experimental) physical 
interpretation of this parameter remains open. 



4-2. Range of validity of GC symmetry 

Armed with the detailed knowledge of the phase diagram for our single-site model, we 
now return to the original question of the validity of the fiuctuation theorem. For the 
GC symmetry to hold for currents of magnitude j, we require that both j and —j are in 
phases B or C (in which e(A) is given by the lowest eigenvalue of H). This immediately 
implies that the symmetry is seen only in the restricted interval [— jmax, jmax] where 

jmax = min(jb, -je). (53) 

This is shown as the shaded regime in figure [5l Note that jmax depends on the initial 
distribution and that the GC symmetry is not seen at all above some critical value of x. 

Now let us investigate the ratio of probabilities for forward and backward currents 
outside the symmetry regime. For non-zero x, then one sees a crossover to 

e(-j) - e(j) = /-j(a,7) +/9(1 - x) + 5{1 - x'^) - jinx - /j(a,7) 

= /3(l-a;)+5(l-a;-i)+jln|^-^^ . (54) 

In other words, for large j, the quantity e(—j) —e(j) is linear with slope In ( — ). At first 
glance, this may appear to be the extended fiuctuation theorem of section [2^61 However, 

sections obtained via Legendre transform do yield the correct form for e(j) in the transition regimes 
(i.e., the large deviation function is still convex). For a mathematical account of the subtleties of 
large deviation theory the reader is referred to [4^; the application to statistical mechanics is discussed 
in [mill]. 
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the slope of the linear section is not as predicted there — in particular, it is not zero for 
an initial steady-state distribution (i.e., the ratio of forward and backward currents does 
not saturate to a constant value). We emphasize also that e(j) itself does not become 
linear for large \j\ as it does in the case of boundary terms which are independent of 
the bulk contribution. It would be interesting to develop a general argument to predict 
the behaviour in cases, such as this, where correlations are important. 

4.3. Generalizations 

Although the calculations of appendix [A] and section 14.11 were for a single-site model 
with the rates Wn = 1 and a Boltzmann initial distribution (l43l) we argue here that they 
also have implications for more general single-site models. 

Firstly, we note that results for the case Wn = a where a is any finite positive 
constant are trivially given by rescaling /3 — > a/3 and 7 — > 07. More generally, we argue 
that for the long-time behaviour only the large-n asymptotics of Wn are relevant and 
thus the results for any bounded W7„ function lim„^oo Wn = a are obtained by the same 
rescaling of /3 and 7. To see this note that the lowest eigenvalue in the gapped state (1401) 
is independent of Wn and the condition for occurrence of instantaneous condensates 
(divergence of the scalar products) depends only on the behaviour forn — 00 as does the 
asymptotic current distribution out from such an instantaneous condensate. However 
the convergence to the long-time limiting behaviour is expected to depend on the form 
of Wn making direct comparison of finite-time simulation results difficult. 

By a similar argument, we expect any initial distribution with the same large-n 
behaviour to lead to the same current large deviations. Since the elements of the lowest 
eigenvector lipo) fall off exponentially with n, this means that any initial distribution 
with super-exponential decay should give the same result as the empty initial site (or any 
other fixed initial configuration), i.e., x = 0. Similarly, if the weight of configurations 
in the initial state decays slower than exponentially, then this corresponds to the case 
X = 1. 

Finally, we remark that results for the current across the output bond can always 
be obtained by the replacements: a ^ 5, /5 <-> 7, p ^ g, A ^ —A, j ^ —j, i.e., by 
left-right refiection. In section [6] we discuss extensions to larger systems. 

5. Numerical results for e(A) 

In [1], the analytical prediction for e(j) was checked via direct Monte Carlo simulation. 
However, it is difficult to get high quality data for long-times since one is looking for 
exponentially unlikely events. In this section we demonstrate instead the application of 
a recent algorithm by Giardina, Kurchan and Peliti [2] to calculate e(A) directly. 
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5.1. Cloning algorithm 

The method is based on an alternative interpretation of ( fT2ll . Here H acts as the 
generator of some kind of time evolution. However, it cannot be an ordinary stochastic 
evolution operator since it does not satisfy the normalization condition ([8]), therefore 
it does not conserve probability. The idea is to consider an ensemble of N identical 
systems (clones), where the number of "individuals" being in the same configuration a is 
denoted by P^. Here the size N of the "population" can vary in time, which means that 
there is no conservation law for the vector \P) as it is not a probability vector anymore 
((s|P) 7^ 1). The off-diagonal elements of H give the currents from one configuration 
to another. The fact that H is not normalized means that the sum of the currents 
from a specific configuration into any other is not necessarily equal to the loss in that 
configuration. This means that individuals can reproduce (i.e. introduce another copy 
of the system prepared in the same configuration as the ancestor by increasing by 
one) or die at given rates, where this rate depends on the configuration a, and is given 

by Ea' ^a'a- 

The process described above can easily be simulated by Monte Carlo methods. 
The quantity to measure is the average rate of growth of the whole population for large 
times, which directly gives e(A). In practice however, one chooses a sufficiently large 
which is kept constant by renormalizing the size of the population after each "birth" 
or "death" (while keeping a record of growth). The method was originally introduced 
for discrete-time update (for details see [2]) but is straightforwardly modified to the 
continuous time case (see also |44j). 

Note that H does not give the full probabilistic description of the above-defined 
stochastic cloning process. The deterministic time evolution (often called rate equation) 
given by H refers only to the averages in a large population [N oo). Therefore, in 
order to obtain reliable results one has to reduce the possible fiuctuations of the (cloning) 
process by choosing large N . 

In a modified version of the algorithm the mean (integrated) current is measured 
instead of the growth rate. This gives the derivative of e(A), which can be then 
numerically integrated (with the fixed condition e(0) = 0). We note that, in this 
algorithm, descendants of a given clone inherit not only the configuration of the parent 
but also the actual value of the parent's integrated current. The advantage of this 
modified algorithm is that it gives less noise. 

5.2. Results for the one-site ZRP with empty initial condition 

We performed the above programme for the one-site ZRP. Since e(A) is known exactly, 
this should be considered as a test of the cloning method. Results are shown in figure 
[6](a) for the case of an empty initial site. It can be seen that the resulting data points 
for N = 10^ and t = 10^ lie very close to the exact {t = oo) results in phases A and B. 
However, in phase C there is a significant deviation and the difference decreases with 
increasing A^. A systematic analysis of the A^-dependence of the cloning method is still 
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Figure 6. Cloning simulation results for the one-site ZRP with empty initial condition 
compared to the exact analytical results. The plots show the result of two versions 
of the algorithm, both with two different values of N. The original cloning algorithm 
(a) produces more noise in the gapless regime than a variant of it (b) whereby the 
derivative of the function e(A) is measured and then numerically integrated to get 
e(A). In the gapless case the convergence in N is slow. For reasonable values of N 
there is a systematic deviation (overestimation) from the exact result. Parameters: 
a = 7 = (5 = 0.1, (3 ^ 0.2. 



lacking so it is not obvious how to set the value of N in general to provide a good 
balance between accuracy and simulation speed. 

Our numerical results suggest that in the gapless phase the simulation could be 
very sensitive to the value of N and could deliver unreliable results for the accessible 
range. We believe that this limitation of the cloning method is related to the gapless 
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spectrum and would not show up in cases where the state space is finite. 

Figure [6](b) shows the results of a modified cloning algorithm, where the derivative 
of e(A) was measured directly and then numerically integrated. Despite the smoother 
results the same type of discrepancy shows up in this version of the simulation. 

In the stationary state of the cloning process (if it exists), the distribution of the 
occupation number of clones is expected to follow the ground-state of H. We measured 
this distribution in a representative point of phase B (A = 1.0, x = 0) and C (A = 2.2, 
X = 0). Results are shown in figurellHHl The measured distribution in phase B follows the 
one suggested by the theory. In phase C (gapless phase) however, there is a significant 
difference between the theory and the measurement. This suggests that in the gapless 
phase the algorithm fails to find the true ground-state, hence e(A) is systematically 
overestimated in the measurement. This is consistent with the results shown in figure 
[6l We note that in phases A and D there is no steady state of the cloning process. 

As discussed above, the cloning method strongly relies on the assumption that the 
number of clones is infinite. It is highly non-trivial to determine what type of correction 
appears if N is finite (as in computer simulations). Our numerical results show that this 
correction is much higher in the gapless phase. This is possibly due to the fact that with 
a finite number of clones one cannot recover the exact distribution but some fluctuations 
are introduced. In the case of a gapped spectrum the system has a strong tendency to 
the ground state, therefore these fluctuations are suppressed and do not play a crutial 
role as long as they are sufficiently small (iV is sufficiently large). In the gapless case 
however, the relaxation to the true ground-state becomes slow and therefore, due to the 
fluctuations, states other than the ground state are also represented in the sample with 
a non-negligible weight (see figures [THHI) . 

It is also possible to define e(A) for finite time, which we denote by e{\)t. 

e(A), = -^ln(e-^-^) (55) 

Since measurements can be made only for finite times, the corrections to e(A) are of 
great interest. It is relatively easy to calculate the finite-time corrections in the single- 
site model based on (1A.13I ). In phases A-B-D the leading contribution comes from a 
pole which leads to an 0{l/t) correction in e(A)t. In the gapless phase C the leading 
contribution is given by the saddle point integration. Since the first order term vanishes 
here the dominant contribution gives 0(t-2/3exp(-e(A)t)) in ((ATSD, which leads to 
a (2/3) ln(Ct)/t correction, where C is some constant. We performed exact numeric 
calculations of e(A)t by evaluating the integrals in ( 1A.13I 1 and this is in good agreement 
with our analytical findings. For details see figure [9l 

5.3. Non-empty initial condition 

It is, in principle, possible to start the system from any initial condition in the cloning 
algorithm. Ideally, one would hope that the method reproduces the exact results in this 
regime. However, since the number of clones is finite, the initial distribution has always 
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Figure 7. Distribution of the occupation number in the "steady state" of the cloning 
algorithm at A = 1.0 (phase B) compared to the analytical results (ground-state of H). 
Parameters: a = 7 = (5 = 0.1,/3 = 0.2, t = 10000, N =10^. 




Figure 8. Distribution of the occupation number in the "steady state" of the cloning 
algorithm at A = 2.2 (phase C) compared to the analytical results. Taking the Hmit 
of expression (|A.3P at A: ^ together with l|A.7P one obtains {n\tpk-^o) ^ 0^"(n + 
y/{y—l)). The algorithm seems to fail in this gapless region and the distribution does 
not converge to the ground-state of H{X). Parameters: a = j = S = 0.1, /3 = 0.2, 
t = 10000, N = 10^ 



25 



exact numeric 




Figure 9. Plotted is the finite-time correction Ae(A)t = e(A)t — e(A) against t in the 
gapless phase at A = 2.2. Exact numerical calculations are compared to the simulation 
results. The cloning method (with fixed N) becomes unreliable for large t. Parameters: 
a = 7 = (5 = 0.1, (3 = 0.2. 

a finite (although arbitrarily large) cutoff. Although this is a minor change in the initial 
distribution, it becomes crucially important when measuring large deviations. A long 
time measurement for an initial distribution with a finite cutoff gives the same results 
as for an initially empty lattice. For an initial distribution with exponential tails (the 
case that we consider) one would need ~ exp(t) number of clones in a measurement 
of length t. This is practically unreachable for reasonably large t. For this reason it is 
not surprising that the algorithm breaks down in the initial-state-dependent phase. 

6. Larger systems 

In this section we extend our discussion to larger systems in order to demonstrate 
the generality of the fiuctuation theorem breakdown for the current fiuctuations. We 
consider the L-site zero-range process with parameters as defined in section [3] and take 
again w„ = 1 (but expect qualitatively the same results for any bounded Wn)- For 
definiteness we assume that the boundary parameters have been chosen to give a well- 
defined steady state with mean ("forward") current to the right. Once again, we focus on 
the behaviour of current fiuctuations across the input ("0th") bond but indicate also how 
to extend our approach to treat the fiuctuations across bulk bonds. In section 16.11 we 
first outline our general method, in particular the use of a powerful mapping to effective 
one-site systems. This leads to a statement about the regime of validity of the GC 
symmetry (section [6^21) and some comments on generalizations (section [Oil . Finally, in 
section 16.41 we give explicit results for a two-site system and make comparisons with 
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numerical data. 



6.1. General approach 

We remind the reader that, for a ZRP of arbitrary size, one can rigorously calculate 
the lowest eigenvalue eo in the gapped phase of the modified Hamiltonian measuring 
the input current together with the associated left and right eigenvectors {ipol and \ipo)- 
Details of the calculations are given in [36], pertinent results are summarized in section [3] 
above. In order to determine the limit of the GC symmetry one has to ascertain the 
values of A at which the scalar products (s|?/'o) and (^ol-Po) diverge whilst keeping in 
mind that the divergence of the normalization of the gapped eigenvalue (^ol^o) indicates 
the cross-over to a gapless regime. To determine the behaviour of current fluctuations 
when these scalar products diverge we can appeal to the heuristic argument based on 
instantaneous condensates. 

This process is considerably simplified by relating the current fluctuations in an 
L-site system to those in a one-site system with effective parameters. We note that 
with the definitions 

di = 7 77 — ., , (56) 

(p-g + 7)(p/?) -7 

^ P{p/q)^~\p - q) , . 

^' p-q-[5 + P{p/qY-' ^ ' 

lip - g) (.o\ 

[P- q + Vip/qy -7 



p - q - (3 + (3{p/qy 

the lowest eigenvalue (l40l) in the gapped phase of the modified Hamiltonian for an L-site 
system (l39l) can be written in the form 



Si = Z : a , a^^,.^L-P (59) 



eo = -^(l-e-) + -^(l-e^), (60) 

Pi + li Pi + li 

and the corresponding left and right eigenvectors are defined by fugacities 

aiQ-^ + Si 

= R (^^) 

Pi + ll 

= -jyrv- ^''^ 

In other words, the one-site marginal for site / and the eigenvalue have exactly the same 
form as in the one site problem with effective ai and ji (A and 6i) determined by all the 
rates to the left (right) of/. This fact relies on the product-state character of the lowest 
eigenvector of the modified Hamiltonian corresponding to the current being measured. 
Notice that, by construction, ai = a, 'Jl = 1, Pl = P and 6l = S. We also remark that 
the parameter combinations aiPi/iPi + 'ji) and ■yiSi / {Pi + 'ji) are independent of / and are 
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equal to the left/right hopping rates of an particle bound state in the corresponding 
exclusion process. Equivalently, in current space we have 

«W = //^.^). (63) 

This important observation allows us to utilize the exact results already obtained for 
the one-site case to build up the phase diagram for an L-site system{Jj] 

Recall that A = corresponds to the mean current. The asymptotic probability of 
seeing a current fluctuation larger than the mean is given by the behaviour of e(A) for A 
negative. To determine the A < part of the phase diagram we carry out the following 
procedure: 

(i) Start with A = and e(A) given by cq of (l60l) . 

(ii) Decrease A until {s\i/jq) diverges on one of the sites, which we label h. From the 
one-site picture we immediately see that this will happen at 

Ai(/i) = Xi{ai,,(3i^,-fi^,6iJ (64) 

with Ai defined as in Table [H Physically, we argue (just as in the single-site 
case) that this divergence corresponds to the "piling-up" of particles on site li. For 
A < Ai(/i) the current fluctuations (across the input bond) only depend on the part 
of the system to the left of site h and we thus expect a crossover to 

e(A) = az,(l-e-^) +7^,(1 -e^), (65) 

which corresponds to a current large deviation function. 

e(j) = /i(«ii>7«i)- (66) 
Just as in the single-site case, eo(A) is continuous but not differentiable at Ai(/i) so 
the two phases in current space are separated by a linear transition regime whose 
explicit form is simply obtained from the effective single-site picture. 

(iii) Now, with such an "instantaneous condensate" on site /i, the left-hand part of the 
system looks just like a system of size h — 1 with right-hand boundary rates p and 
q. We can then write down the new ground state \ipo) in terms of redefined effective 
parameters a/, /?/, •ji and 6i. 

(iv) Repeat steps ([n]) and ([ml) recursively until... 

(v) At some value of A, (s|V^o) (with lipo) a function of the effective parameters) diverges 
on site 1 and the large deviation function takes the form 

e(A) = a(l-e-^)+7(l-e^) (67) 
e(j) =/,(a,7)- (68) 
(again with a linear transition regime in j-space). In other words, for very large 
forward currents, only the hopping parameters across the input bond are significant. 

II Although our results can be applied to large systems by taking the thermodynamic limit L oo, 
this limit does not necessarily commute with the long-time limit t ^ oo. Our analysis therefore does 
not address the form of current fluctuations in a genuinely infinite system (taking the limit L ^ oo 
first, followed by t ^ oo). 
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In a similar fashion one can investigate what happens for currents smaller than 
the mean by increasing A from zero. For simplicity, let us first take a fixed initial 
configuration so that {iPq\Po) is always finite. At A2(/) = A2(a;, ji, 5/), (ipoli^o) diverges 
on site /, corresponding to a transition to a gapless phase. The minimum value of A2(/) 
occurs at some I2 (which depends on the parameters of the model) so we first see a 
crossover to 

e(A) = ai, + + Ih + K - 2\/(a«2e"^ + (5/J(A2 + Ih^^) (69) 

e(j) = fjiah^lh) + fjif^h^^h)- (70) 

In other words, to sustain a large backwards current we have an instantaneous 
condensate on site I2 and the product of current distributions across the two independent 
parts of the system. For increasing backwards currents, one can then repeat the 
procedure for the two subsystems (each with redefined effective parameters) , increasing A 
and looking for the next site where (^ol^o) diverges. Eventually, for very large backward 
currents we expect instantaneous condensates on all sites and a current large deviation 
function given by 

e(j) = /,(«, 7) + (^ - 1)/,(P, q) + fjiP, S) (71) 

For a distribution over initial configurations, the situation is more complicated but 
again we can make progress based on the effective one-site picture. For an initial particle 
distribution which is a product measure of Boltzmann distributions with site dependent 
fugacity xi, then {i/jqIPq) will diverge on the Ith site at Xa^I) = \4{ai, Pi,'yi,6i,xi) 
leading to a crossover to an initial-state dependent regime. If I4 is the (initial-condition- 
dependent) value of / corresponding to the minimum of A4(/) then, for large initial 
fugacities we expect to find first a transition to 

e(A) = ai, + (3i, + -fi, + 61, - (A4 + - {ai^e'^ + 6u)xl^^ (72) 

e(j) = /j(az4,7u) + A4(l -a^u) + 5u(l -x^^) + jlnx;4. (73) 

Once again, we see a factorization of the asymptotic current distribution across the 
two sub-systems to the left and right of ^4. In principal, one can try to build up the 
complete phase diagram by next checking for the divergence of {%Ijo\Po) on other sites 
as A is further increased. Another scenario, for small initial fugacities, involves first a 
transition to a gapless phase and then a transition to the initial condition dependent 
phase when A = A3(/) as a function of the redefined effective rates. 

6.2. Range of validity of GC symmetry 

The effective one-site picture provides an elegant way to summarize the range of validity 
for the GC symmetry in ZRPs of arbitrary size. By comparison with (i53ll and with the 
shaded regime in figure [5l we see that, for given initial fugacities {x}, the symmetry 
relation for input current will only be obeyed in the following current range: 

\i\ < mm{jb{ai^,(3i^,'yi^,6i,),-jeia^,l3i^,'yu,6^,xij) (74) 
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where /« is the value of / at which |Aj(/)| takes its minimum value as a function of the 
effective rates given in ([56l) - (l59|l . In the case where je is positive for some / then the 
symmetry will not be observed at all. 



6. 3. Further generalizations 

We remark that the same condition as ( l74l l should also apply for ZRPs with quenched 
disorder (i.e., with p and q bond dependent) — one needs only to determine the 
appropriate form for the effective rates. Indeed, building on the recent determination 
of the stationary state for such systems [45], we can show that the relevant parameter 
combinations are 

nl-l 2l 
^4=0 llfc=l qi_^ 

_ 1 

~ ^i-i 1 ffi pTT v ' / 

^i=0 llfc=l g^^^ 

ttL — i Ql+i 

where pi (qi) is the right (left) hopping rate across bond / and by definition po = a, 
Qo = It Pl = P and = 6. One can state (I75II78P alternatively as a simple recursion 
relation: 

= (80) 

S,_, = (82) 

A + qi-i 

Not surprisingly, these are again the renormalized two-particle hopping rates. 

We conclude this subsection by indicating how to extend our approach to treat 
the current fluctuations across other bonds. In this case, the effective parameters are 
unchanged, but in all sites to the left of the bond in question the fugacities (l6Tll and (l62l) 
must be replaced by 

= ir^ («^' 

Pi + li 

h = ^-^P^. (84) 
A + 7/ 

corresponding to the expressions for the output bond of a single-site system. 
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6.4- Example: Two-site system 



In principle, the general approach outlined above can be used to determine the behaviour 
of current fluctuations in ZRPs of arbitrary size. However the implementation rapidly 
becomes tedious and the phase diagram complicated. Here, as a simple test case, we 
present results for a two-site system with empty initial condition (xi = for all /). 

It is obvious that for any choice of parameters obeying a — •yKp — q</3 — 6, the 
system has a well-defined steady state (i.e., no boundary condensation). One also gets 
a stationary state if the conditions 

a-^<(3-5<p-q (85) 

and 

ap jq 



</?-7 



(86) 



p + 7 P + 7 

are obeyed. These inequalities can easily be understood in the exclusion picture where 
the first two particles form a bound state with hopping rates and The condition 
for the bound state of this and the remaining single particle is just (l86ll . 

For definiteness, in the remainder of the discussion we assume that the rates 
obey (l85l) and (l86ll . In this case, as A is decreased from 0, {slipo) diverges first on 
site 2 (i.e., h = 2). Similarly, as A is increased from 0, {ipoli^o) also diverges first on site 
2 (i.e., I2 = 2). Appealing to the heuristic "instantaneous condensates" picture we then 
obtain the following different regimes for the current large deviation function. 



e(j) 



/i(«,7) 

P+7/ 

A 7g+/3p+/37 ' 'yq+(3p+f3j J 

/,(a,7) + /,(p,g) + /,(/3,5) 



jy&q__\ 



A2 
Al 
B 

CI 
C2 



(87) 



with 



A2 
Al 
B 
CI 
C2 



jaia,p,-f,q) < j 

3. 



( "P 



Vp+7' ' P+7 

Vp+7' ' ' p+7' 
3c{a,p,-f,q) <J <3c(:^,P,f^ 



<3<3b[^,P, 



6 



^,6 

7' 



(88) 



3 > Jc(a,P,7,g) 

Here the currents ja, jb and jc are as defined in Table [1} we label the phases in analogy 
to A, B and C in the single-site case and note that there are intermediate transition 
regions at the A2-A1 and Al-B crossovers. From the argument in 16.21 it is clear that 
the Gallavotti-Cohen symmetry should hold for small currents j < jb (^^) f^i ^5 ) • 
Finally, we compare the Legendre transform of (l87l) with results for e(A) obtained 
via the cloning algorithm (see section [5|). As shown in figure [TOl we find excellent 
agreement except for A > A2 (^^! ^) where e(A) is given by the lowest eigenvalue 

of a gapless spectrum. Just as in the single-site case we attribute this to limitations 
imposed by the finite number of clones. 
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Figure 10. Cloning simulation results for two-site ZRP with empty initial condition 
compared with exact analytical results. Agreement is excellent except in the gapless 
phase. Parameters: a = j = 6 = q = 0.1, /? = 0.15, p = 0.24, t = 5000, N = 10". 



7. Summary 

This paper focuses on the Gallavotti-Cohen (GC) fluctuation symmetry in the context 
of continuous-time Markov processes. In particular, it offers a contribution towards 
understanding the potential breakdown of this symmetry in systems with infinite state 
space. To set the scene for this, we flrst discussed how the symmetry is manifested for 
the large deviation function of currents in systems with finite state space. Our proof 
was based on that of Lebowitz and Spohn [H] (or Derrida et al. [46j) but slightly 
extends their work by considering more general currents. In the notational framework 
of this paper, the action functional of Lebowitz and Spohn can be obtained by choosing 
Qa,a' = ln(W(j',cr/W(j,cr') and the GC symmetry then holds for J with E = 1. In this case, 
as shown in pTj, J can be interpreted as an entropy current. 

The above-mentionned proof of the fluctuation relation holds only for systems with 
finite state space — in systems with infinite state space the GC symmetry can be broker^. 
In the spirit of van Zon and Cohen's earlier work for Langevin systems [23l [17], we 
next recapitulated the ideas leading to an "extended fluctuation theorem" in which 
(for stationary initial state) the ratio of probabilities of given forward (j) and backward 
(— j) currents approaches a constant for large values of j. We argued that this particular 
form of breakdown relies on the independence of bulk and boundary contributions to 
the current. The central aim of this paper was to study analytically a simple model 
where this condition is not met — the partially asymmetric Zero Range Process (ZRP) 

% We remark that an infinite state space is not a sufficient condition for a breakdown of the Fluctuation 
Theorem. A counterexample is the ZRP with unbounded w„ as explained in section [321 
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on an open lattice [35l 136] . 

Specifically, we expanded on earlier results in [1], and showed in detail how to 
calculate the large deviations of the average particle current j in the one-site ZRP. 
The phase behaviour of the large deviation function can be physically understood by 
considering the possibility for an arbitrarily large number of particles to pile-up on the 
site ("instantaneous condensation"). Even in this single-site case, the phase diagram 
(parametrized by the time-averaged current and the initial state) shows a complex 
picture with first and second order phase transitions. These transitions are analogous 
to ordinary equilibrium phase transitions since formally the current large deviation 
function can be considered as a kind of dynamical free energy functional. In general, 
the GC symmetry holds only in a restricted interval [— jmax, Jmax], where jmax depends 
on the initial state. The relevant parameter of the initial state, which is denoted by 
X, can be considered as a kind of initial temperature. For large values of x the regime 
of validity of the fiuctuation relation shrinks (jmax decreases) and reaches zero at a 
finite X. According to the above interpretation, this means that above a critical initial 
temperature the fiuctuation theorem does not hold, even for small currents. 

Although the one-site ZRP might be thought to be oversimplified, it is very 
instructive and already exhibits most features of the multiple-site version, which we 
have also studied in detail. In particular, we developed a heuristic argument based on 
mapping to an effective one-site system and considering the appearance of instantaneous 
condensates. The phase diagram can, in principle, be derived for any number of sites 
but it becomes increasingly complicated. As an example we showed how the current 
large deviation function can be obtained for a two-site ZRP. Again the phenomenon of 
instantaneous condensation is the physical mechanism leading to a breakdown of the 
GC fiuctuation theorem. 

Outside the symmetry regime, we found that the behaviour is somewhat different to 
that predicted by the extended fiuctuation theorem of van Zon and Cohen, presumably 
due to the correlations (between bulk and boundary terms) in our model. There are 
evidently still open questions relating to the characterization of the fiuctuation behaviour 
of different systems beyond the limits of Gallavotti-Cohen symmetry. For example, in 
another recent work [47j it was shown that the distribution of work fiuctuations for 
a Brownian particle with Levy noise (i.e., infinite variance) do not even have a large 
deviation form. It might be interesting to look for analogues of the resulting "anomalous 
fiuctuation theorem" in the present framework of many-particle dynamics described by 
a master equation. 

Our analytical results also made it possible to test the recently proposed cloning 
algorithm [2] that measures the Legendre transform of the large deviation function 
directly. The strength of this method is its efficiency and relative simplicity. The 
algorithm reproduces our analytical results in several phases. However, there are notable 
deviations in a phase where the spectrum of H, the effective Hamiltonian governing the 
current fiuctuations, becomes gapless. Another weakness of the algorithm is that it is 
unable to capture the initial-state-dependence of the fiuctuations. It is only the case 
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of a fixed initial configuration that this method can safely be used. These observations 
call for the development of a novel numerical algorithm which would not break down 
in cases where the cloning method does. One possibility here would be to develop a 
transition path sampling method. 
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Appendix 

A. Analytical calculations for single site model 
A.l. Spectrum 

Naively, one expects that the large deviation of current fluctuations is given by the 
lowest eigenvalue of H according to (fT2l) . For this reason, we here study the spectrum 
of (l44l) . First we transform H into the symmetric form 



H'{X) = <^H{X)<I>-^ 

using the diagonal operator 

f (p 




(A.l) 



with 



ae^^ + 5' 



(A.2) 



0^ 

V : : : • 

The diagonalization of H'{X) is easily done by a Fourier transformation. We introduce 



y{X) = \J + (3) (ae ^ + 5)/ (/3 + 7) which determines the character of the spectrum 
as follows: 

y{X) > 1 Here the spectrum is entirely continuous. The eigenvector |^'(A;)) 
corresponding to wave number k G (0, vr] takes the form 

\^'{k)) = \l=-Y,sm{kn + <f)\n), with e'^^ = ' 



n=0 



l/(A)e- 



-ik 



(A.3) 



and the corresponding eigenvalue is 



e{k) = a + /? + 7 + 5- 2 J(7e-^ + (3) (ae"^ + 6) cos k. 



(A.4) 
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y{X) < 1 In addition to the above continuous spectrum a discrete band appears with 



|^'(0)) = Jl-y^E?/"!^) ande(0)=a + 5-(/3 + 7)l/'. (A.5) 



n=0 



Here this is the lowest eigenvalue, correspondingly the spectrum becomes 
gapped. Notice that e{k ^ 0) 7^ e(0). 

We note that it is easy to prove that the above set is complete, i.e., 

^n,-{m\ij'{k)){^'{k)\n)dk y{X)>l 

\ Jo{m\i;'{k)){tP'{k)\n)dk+ {m\^'iO)){ij'iO)\n) y(A) < 1 ^ ' 

The right and left eigenvectors and of H{\) can be obtained by 

applying the operator $: 

for k G [0, tt] . In summary, for the lowest eigenvalue eo(A) (infimum of the spectrum) 
of H{X) we obtain 

c^ + S-i(3 + i)y{Xf y{X)>l 
°^ ^ 1 a + /3 + 7 + 5-2(/3 + 7)2/(A) y(A) < 1 ' ^ 



It can easily be seen that this expression satisfies the symmetry relation (l2Tll . since 
y{X) = y{E - A) with = (a/3) / (7^). 

A. 2. Calculation of the large deviation function 
As a first step we write 

{s|e-^*|Po) = (l-x)5]5]x"(m|e-^», (A.9) 

m n 

where {m\ is a row vector with a T at the mth position and 'O's elsewhere. Inserting a 
complete set of eigenvectors we can write this in the form 

(s|e-^*|Po) = (1 -x)EEx" r {m\i,{k)){iP{k)\n)e-'^^^' dk 

rn n 

+e (1 - y) (1 - x) E Ex'^(m|^(0))(^/'(0)|n)e-^(°)* (A.IO) 

m n 

Here 6 denotes the Heaviside function. Using the k —k symmetry of the eigenvectors 
and eigenvalues in the contribution of the continuous part, the above integral can be 
rewritten in the form 

JO 27r Jo ^ ^ 

After the substitution z = and using (IA.3I 1 we obtain 

(s|e-^*|Po) = (1 - a;) E E / Ln-m-l ^ V^^^n+nA ^-e{z)t 

m n J|2|=l V Z — y J 

+e (1 - y) (1 - x) E E (1 - 2/') j/"+™e-^(°)*, (A.12) 
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where by e{z) we mean e{k{z)) based on the above substitution. In order to be able 
to perform the infinite sums in (IA.12I) . we have to choose the contour of the integral 
carefully. In the first term of the integral there is a pole only at z = 0. Here we chose 
the contour to be a circle of radius < \z\ < (0a;)~^ (denoted by Ci). In the second 
term we deform the contour to run along a circle around the origin with infinitesimal 
radius (denoted by C2). Notice that doing so, for ?/ < 1 we pick up a pole contribution 
at z = y, which just cancels with the last term in (IA.12I) (since e(0) = e{y)). Finally we 
obtain 

2.ix0/c.(.-^)(.-i) 
1 — X f {yz — 1) e 



(s e 



, , , -dz. (A.13) 

2^1^ Jc. {z -^^){z-<p){z-y) 

To obtain the large deviation function we need to study the limit of this integral as 
t — s> 00. This enables us to build up the phase diagram as shown in section [4Tl 
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